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Abstract 

Modifications to an existing three-dimensional, im- 
plicit, upwind Euler/Navier-Stokes code (CFL3D Version 
2.1) for the aeroelastic analysis of wings are described. 
These modifications, which were previously added to CFL3D 
Version 1.0, include the incorporation of a deforming mesh 
algorithm and the addition of the structural equations of mo- 
tion for their simultaneous time-integration with the govern- 
ing flow equations. The paper gives a brief description of 
these modifications and presents unsteady calculations which 
check the modifications to the code. Euler flutter results 
for an isolated 45° swept-back wing are compared with ex- 
perimental data for seven freestream Mach numbers which 
define the flutter boundary over a range of Mach number 
from 0.499 to 1.14. These comparisons show good agree- 
ment in flutter characteristics for freestream Mach numbers 
below unity. For freestream Mach numbers above unity, the 
computed aeroelastic results predict a premature rise in the 
flutter boundary as compared with the experimental bound- 
ary. Steady and unsteady contours of surface Mach number 
and pressure are included to illustrate the basic flow charac- 
teristics of the time-marching flutter calculations and to aid 
in identifying possible causes for the premature rise in the 
computational flutter boundary. 

Nomenclature 

Aij generalized aerodynamic force resulting from 

pressure induced by mode j acting through 
mode i 

b root semichord 

c root chord 

Ci generalized damping of mode i 

C p pressure coefficient 

k reduced frequency, 

k{ generalized stiffness of mode i 
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m measured wing panel mass 

rrii generalized mass of mode i 

Moo freestream Mach number 

q { generalized displacement of mode i 

Qi generalized force of mode i 

T dimensional time 

Ui load vector 

U F flutter speed 

Uqo streamwise freestream speed 

v volume of a truncated cone having 

streamwise root chord as lower base 
diameter, streamwise tip as upper base 
diameter, and panel span as height 

Xi state vector 

a steady-state angle of attack 

0 integral of the state-transition matrix 

M mass ratio, 

<jj damping associated with mode j of the 

aeroelastic response 

$ state-transition matrix 

w angular frequency 

ufj frequency associated with mode j of the 

aeroelastic response 

u> a uncoupled natural frequency of the wing first 

torsion mode 

Introduction 

Research during the last decade on the application of 
computational fluid dynamics (CFD) methods to unsteady 
flows and aeroelastic analysis has been rapidly progressing. 
Edwards and Malone 1 recently presented a survey on the 
status of computational methods for unsteady aerodynamic 
and aeroelastic analysis with an emphasis on methods for 
transonic flows. The transonic speed range has been a maijp 
focus of activity because flutter dynamic pressures re tjjp^ 


ically critical (lower) in this speed range. Much of this re- 
search, especially for three-dimensional configurations, has 
focused on the development of finite-difference methods for 
the solution of the transonic small disturbance (TSD) and 
full potential (FP) equations. 1 One reason for the focus 
on the FP and TSD methods is that the reduced memory 
and run-time requirements of these methods in comparison 
with the higher-order methods have made them more viable 
for use in aeroelastic analyses of three-dimensional config- 
urations. Edwards and Malone 1 reported on 13 aeroelastic 
studies of flexible wings and flexible wings/rigid body con- 
figurations which utilized the TSD and FP methods. These 
studies, which compare flutter boundary calculations with 
experimental data, provide important applications of these 
CFD methods. 

Reference 1 points out that many critical challenges 
facing computational aeroelaticity will require the modeling 
of increasingly more complex flow physics. To meet these 
challenges, researchers have begun to develop higher-order 
methods involving the Euler and Navier/Stokes equations for 
unsteady aerodynamic and aeroelastic analysis. With recent 
advances in algorithm development and computer hardware, 
higher-order methods utilizing the Euler and Navier/Stokes 
equations have been used in three-dimensional aeroelastic 
applications; 2-11 however, the number of these applications 
lags behind those utilizing the TSD and FP methods largely 
in part because of their increased computational require- 
ments. 

The research described in Refs. 2-11 represents impor- 
tant steps in the development of three-dimensional Euler and 
Navier/Stokes methods for aeroelastic analysis. However, 
continuing studies are needed to validate these methods for 
the prediction of aeroelastic response and flutter. In Ref. 3, 
Robinson et al. performed time-marching flutter calculations 
for an isolated 45° swept-back wing using an Euler code, 
A novel aspect of the capability described in Ref. 3 was 
the deforming mesh algorithm which was used to move the 
mesh so that it conformed continuously to the instantaneous 
position of the wing. The results presented compared fa- 
vorably with the experimental data and with results from a 
transonic small disturbance code for the single flutter point 
analyzed. The purpose of the present work is to further 
demonstrate and assess the capability presented in Ref. 3 by 
completing the flutter boundary for the simple, well-defined, 
isolated 45° swept-back wing configuration using the Euler 
equations. The wing analyzed in these studies is the first 
AGARD standard aeroelastic configuration for dynamic re- 
sponse, and its flutter data is an accepted set with which to 
test codes. 12 In Ref. 3, modifications were made to an ex- 
isting three-dimensional, unsteady Euler/Navier-Stokes code 
(CFL3D Version 1.0) for the aeroelastic analysis of wings. 
These modifications included the incorporation of a deform- 
ing mesh algorithm and the addition of the structural equa- 
tions of motion for their simultaneous time-integration with 
the governing flow equations. The deforming mesh algo- 


rithm and the structural equations of motion described in 
Ref. 3 have been added to the most recently released ver- 
sion of CFL3D, Version 2.1, by the present authors. This 
paper gives a brief description of these modifications and 
presents unsteady calculations which check the modifica- 
tions to the code. Results from calculations performed for a 
rigid wing undergoing forced pitching and plunging motions 
are presented to test the performance of the deforming mesh 
algorithm. Aeroelastic results for the 45° swept-back wing 
at a freestream Mach number of 0.9 are compared to those 
presented in Ref. 3 to check the addition of the structural 
equations of motion. Calculated flutter results for the same 
45° swept-back wing are compared with the experimental 
data for seven freestream Mach numbers which define the 
flutter boundary over a range of Mach number from 0.499 to 
1.14. Steady-state Mach contours of the initial flowfields are 
also included in the discussion of the aeroelastic results to 
illustrate the basic flow characteristics of the time-marching 
flutter calculations at selected freestream Mach numbers. In- 
stantaneous surface pressure contours during an aeroelastic 
transient at Moo = 0.99 are presented to demonstrate changes 
in the flowfield which are induced by the aeroelastic motions. 

Upwind Euler/Navier-Stokes Algorithm 


The time-dependent Euler equations are solved within 
the CFL3D 13, 14 code by a three-factor, implicit, finite- 
volume algorithm based on upwind-biased spatial differenc- 
ing. The algorithm, which is a cell-centered scheme, uses 
upwind differencing based on either flux-vector splitting or 
flux-difference splitting. Both types of upwind differencing 
account for the local wave-propagation characteristics of the 
flow and sharply capture shock waves. Also, because these 
schemes are naturally dissipative, additional artificial dissi- 
pation terms are not necessary. Several types of flux-limiting 
are available within the code to prevent oscillations in the 
solution near shock waves. These oscillations are typically 
found in results from higher-order schemes. For unsteady 
cases, the original algorithm contains the necessary met- 
ric terms for a rigidly translating and rotating mesh which 
moves without deforming. For cases involving a deforming 
mesh, however, an additional term accounting for the change 
in cell volume must be included in the time-discretization of 
the governing equations. This modification is implemented 
as described in Ref. 3. The aeroelastic equations of motion 
were implemented in the more recent Version 2.1 because, 
in addition to other improvements, this code contains several 
options for computing multi-block solutions which will be 
utilized in future computations. 

Deforming Mesh Algorithm 


In the time-marching aeroelastic calculations, the mesh 
must be updated at every time level so that it conforms to 
the aeroelastically deformed shape of the wing. Because the 


2 



aeroelastic motion of the wing may be general in nature and 
are not known a priori, a general mesh updating procedure is 
necessary. One such method, the deforming mesh algorithm, 
models the mesh as a network of springs and solves the 
static equilibrium equations for this network to determine 
the new locations of the mesh grid points. This algorithm 
was originally developed by Batina 15 for tetrahedral cells 
and extended by Robinson et al. 3 for hexahedral cells. The 
edge of each hexahedral cell is modeled as a spring whose 
stiffness is inversely proportional to the length of the edge 
raised to a power. In order to control cell shearing and to 
prevent the collapse of the cell, diagonal springs are added 
along the faces of each cell. Similarly, the stiffness of the 
diagonal springs is also inversely proportional to the length 
of the diagonal raised to a power. As suggested in Ref. 3, a 
power of three was used in the present calculations. 

At each time level, the grid points on the outer boundary 
are held fixed, and the displacement of the wing surface is 
specified. For aeroelastic calculations, the displacement is 
determined from the integration of the structural equations 
of motions. The new locations of the interior grid points are 
then determined by solving the static equilibrium equations 
which result from a summation of forces in the x, y and 
z coordinate directions at each grid point. These static 
equilibrium equations are solved using a predictor-corrector 
method. The new grid point locations are first predicted 
by an extrapolation from the previous two time levels and 
then corrected by using several Jacobi iterations of the static 
equilibrium equations. In the present calculations, four 
Jacobi iterations are sufficient to move the mesh. 

Although the deforming mesh algorithm is a general 
procedure, the current implementation in CFL3D is restricted 
to meshes of C-H topology. This limitation is due to the 
specialized treatment required for the mesh boundaries. For 
example, in a C-H mesh the plane of points represented 
by the maximum i index would be an outer boundary; 
however, in a C-O mesh the plane of points represented 
by the maximum i index would be an internal cut. 

Time-Marching Aeroelastic Analysis 


In a time-marching aeroelastic analysis, the calculation 
of each flutter point is begun by obtaining a static aeroelastic 
solution about the wing. There are several ways to obtain a 
static aeroelastic solution. One way is to compute a steady- 
state solution about the rigid wing. The next step is to 
allow the rigid wing to deform to the steady loads until 
a static aeroelastic solution is obtained. The calculation 
of the static deformation is obtained using the aeroelastic 
equations of motion. For this time-marching calculation, 
however, the structural damping of the wing is set to a 
number, around 0.99, such that the dynamic system is near 
critically damped. The aeroelastic equations of motion are 
then marched simultaneously in time with the governing 
flow equations until the wing no longer deforms under the 


aerodynamic loads. Since the system is very damped, this 
calculation does not require a great amount of computational 
time. The static aeroelastic solution is then used as the 
starting point for the time-marching dynamic aeroelastic 
solution. Since the wing analyzed in the present work has 
a symmetric airfoil section, at zero degree angle of attack 
this configuration will have no static deflection. Therefore, 
for this wing, the steady, rigid solutions could be used 
as the starting solutions for the aeroelastic time-marching 
calculations. 

In order to bracket the flutter point, the static aeroe- 
lastic and dynamic aeroelastic computations are computed 
at several values of dynamic pressure (typically three val- 
ues) which ranged from 80% to 120% of the experimental 
values depending on the freestream Mach number. In each 
of the dynamic aeroelastic calculations, the motion of the 
wing is initiated by specifying a small initial velocity for 
the first two modes. The resulting transients are analyzed 
with a modal identification technique for their damping and 
frequency content. The computed flutter dynamic pressure 
and frequency can then be determined by interpolating the 
specified dynamic pressures and the computed frequencies 
to the zero damping value of the dominant mode at flutter. 
The subsequent sections contain a brief description of the 
aeroelastic equations of motion, the time-marching solution 
procedure, and the modal identification technique. 

Aeroelastic Equations of Motion 

The aeroelastic equations of motion that are incorpo- 
rated within CFL3D are similar to those described in Refs. 3 
and 16. In this formulation, the equations are derived by as- 
suming that the general motion of the wing can be described 
by a separation of time and space variables in a finite modal 
series. This modal series consists of the summation of the 
free vibration modes weighted by the generalized displace- 
ments. After applying Lagrange’s equations to this system, 
the aeroelastic equations of motion can then be written for 
each vibration mode i as 

mqi + Ci4i + = Q% 0 ) 

where <?, is the generalized displacement, m % is the general- 
ized mass, Ci is the generalized damping, k{ is the general- 
ized stiffness, and Qi is the generalized aerodynamic force 
computed by integrating the pressure weighted by the mode 
shapes. The superscript dots in Eq. (1) represent differenti- 
ation with respect to time. 

Time-Marching Solution 

The solution procedure implemented in CFL3D for 
integrating Eq. (1) is that described by Edwards et al. 17, 18 
The linear state equations are written as 

Xi = Ax{ + Bui ( 2 ) 
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where A and B are coefficient matrices that result from the 
change of variables Xi = [<? 3 qi] and u { is the nondimen- 
sional representation of the generalized force Qi . Equation 
(2) is integrated in time using the modified state-transition 
matrix structural integrator 17 implemented as a predictor- 
corrector procedure. The prediction for x" +1 , £t n+1 , is 
given by 

£,- n+1 = + QB( 3< - u?- l )f2 (3) 

where 3> is the state-transition matrix and 0 is the integral 
of the state-transition matrix from time step n to n + 1. 
The predicted value of the generalized displacement x; n+1 
is used to update the mesh for the next flow field calcula- 
tion which is used in turn to evaluate the nondimensional 
generalized force u, n+1 . These values are then used in the 
corrector step to determine x? +i , given by 

x* +1 = + 0£(u,- n+1 + u?)/2 (4) 


Modal Identification Technique 

Damping and frequency characteristics of the aeroe- 
lastic responses are estimated from the response curves by 
using the modal identification technique of Bennett and 
Desmarais. 19 The modal estimates are determined by a least 
squares curve fit of the responses of the form 

m 

qi{T) = «o + ^ e° lT [ctj cos(u>jT) + bjsin{uj } T)) 

j = l 

f = 1,2,3,... 

where qi is the generalized displacement of the natural 
vibration mode i (as previously defined) and where <x ; and ujj 
are the damping and frequency, respectively, associated with 
mode j of the aeroelastic response. The number of modes m 
determined in the curve fit of the response is usually greater 
than or equal to the number of modes initially excited. 

Pulse Transfer-Function Analysis 


Generalized aerodynamic forces (GAF’s) can be ob- 
tained by calculating several cycles of a forced harmonic 
oscillation and using the last cycle of oscillation to deter- 
mine the load. This requires one time-marching calculation 
for each value of reduced frequency and each mode of in- 
terest. In contrast, the GAF’s may be determined for a wide 
range of reduced frequency in a single time-marching calcu- 
lation for each mode using the pulse transfer-function anal- 
ysis. In the pulse analysis, the unsteady force is computed 
indirectly from the response of the flowfield due to a wing 
motion that is represented by a smoothly varying, exponen- 
tialFy shaped pulse. A fast Fourier transform of the unsteady 
force is divided by the Fourier transform of the displacement 


to obtain the GAF. The pulse transfer-function analysis has 
been previously employed to determine the GAF’s which 
are used in aeroelastic analyses. 3, 16, 20-22 Results presented 
in Refs. 3, 16 and 20-22 have shown that the analysis is 
valid for predicting the small perturbation response about a 
nonlinear flowfield. 

Wind Tunnel Model Description 


The wing being analyzed in this study is the first 
AGARD standard aeroelastic configuration for dynamic re- 
sponse, Wing 445.6, 12 which was tested in the Transonic Dy- 
namics Tunnel (TDT) at NASA Langley Research Center. 23 
The Wing 445.6 has a quarter-chord sweep angle of 45°, 
a panel aspect ratio of 1.65, a taper ratio of 0.66, and a 
NACA 65A004 airfoil section. A planform view of this 
wing is shown in Fig. 1. Several different models of the 
Wing 445.6 were tested in the TDT including both full span 
and semi-span models. The model used in this study was one 
of the semi-span wind tunnel-wall-mounted models which 
was constructed of laminated mahogany. In order to obtain 
flutter data for a wide range of Mach number and density 
conditions in the TDT, holes were drilled through several of 
the mahogany wings to reduce their stiffness. The aerody- 
namic shape of the original wing was preserved by filling the 
holes with rigid foam plastic. A photograph of a weakened 
model mounted in the TDT is shown in Fig. 2. The model 
designated as “WEAK3” in Ref. 23 is analyzed herein. The 
flutter data for this model tested in air is reported in Ref. 23 
over a range of Mach number from 0.499 to 1.141, 

The Wing 445.6 is modeled structurally using the 
first four natural vibration modes which are illustrated in 
Figs. 3(a) and (b). Figure 3(a) shows oblique projections of 
the natural modes while Fig. 3(b) shows the corresponding 
deflection contours. These modes which are numbered 1 
through 4 represent first bending, first torsion, second bend- 
ing, and second torsion, respectively, as calculated by a fi- 
nite element analysis. 23 The modes have natural frequencies 
which range from 9.6 Hz for the first bending mode to 91.54 
Hz for the second torsion mode. As suggested in Ref. 23, 
the experimentally determined modal frequencies were used 
in the time-marching flutter analysis. For the cases consid- 
ered in this study, no structural damping is included in the 
aeroelastic equations of motion (c; = 0 for all modes). 

Results and Discussion 

Results are presented in this section for calculations 
about the Wing 445.6. All of the computational results were 
obtained using a 193 x 33 x 41 C-H-type grid with 193 
points wrapped around the wing and its wake (129 points on 
the wing surface), 41 points distributed from the wing root 
to the spanwise boundary (25 points on the wing surface), 
and 33 points distributed radially from the wing surface to 
the outer boundary. This mesh topology was chosen rather 
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Figure i Planform view of Wing 445.6. 




Figure 2 Wing 445.6 model in the NASA 
Langley Transonic Dynamics Tunnel. 

than a C-O-type topology because the wind tunnel model 
has a sheared-off tip. This can be seen in the photograph of 
the model in Fig. 2. A partial view of the surface mesh on 
the wing and symmetry plane is shown in Fig. 4. Its outer 
boundaries extend 10 local chord lengths to the upstream and 
downstream boundaries, 10 local chord lengths to the upper 
and lower boundaries, and 1 semi-span length off of the tip. 
Note that this grid is identical to the one used in Ref. 3. 

For all of the calculations, the Euler equations are 
solved using flux-vector splitting and a smooth flux-limiter. 
Convergence to steady-state was accelerated using local 
time-stepping, mesh sequencing and multi-grid cycling. For 
time-marching calculations, the nondimensional global time- 
step (based on the root chord and the freestream speed of 
sound) was 0.05456. The in-core memory requirement for 
this grid was 25 Mw on a Cray-2 computer which is 15.6 
Mw over what is required for the unmodified CFL3D Ver- 
sion 2.1. Recent coding modifications have reduced the ad- 
ditional memory requirements for the aeroelastic version by 
approximately 50%, and so the current grid would now re- 




Mode 3, 48.35 Hi Mode 4, 91 .54 Hi 


(a) oblique projections. 





(b) deflection contours. 

Figure 3 Natural vibrations modes for Wing 445.6. 

quire 17.4 Mw. Each of the steady-state calculations re- 
quired approximately 3-4 hours of CPU time on a Cray-2 
computer to converge the solution to an acceptable level (6 
orders of magnitude), Aeroelastic transients were computed 
at each dynamic pressure for approximately 2 cycles of the 
lowest frequency modal motion. These calculations typi- 
cally required 8 hours of CPU time. Since aeroelastic tran- 
sients were computed for three different values of dynamic 
pressure at each freestream Mach number, the total compu- 
tational cost for each flutter point (including the steady-state 
solution) was approximately 28 hours of CPU time. 

Pulse Transfer-Function Results 

The generalized forces for the Wing 445.6 at M^ = 0.9 
were computed using the pulse transfer-function analysis 
method described in a preceding section. The pulse cal- 
culations were restarted from a steady-state flow condition 
at an angle of attack of a = 0°. A plunging motion and a 
pitching motion about the root quarter chord, which are de- 
fined as modes h and 9 respectively, were analyzed. These 
simple “modes” were chosen in order that the motion of the 
wing could be simulated not only by the deforming mesh 
algorithm but also by a rigid translation and rotation of the 
grid. The maximum amplitude of the plunging motion was 
0.01 root chord lengths, and the maximum pitch amplitude 
was 1°. The results of the pulse analyses, shown in Fig. 5, 
are plotted in terms of the real and imaginary components of 
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Figure 4 Partial view of the 193 x 33 x 41 computational 
grid on the wing surface and symmetry plane. 


the unsteady forces as a function of the reduced frequency k 
which is defined by The generalized force Ahh is the 
lift due to plunge, Am is the lift coefficient due to pitch- 
ing, A$h is the pitching moment due to plunge, and Aee is 
the pitching moment due to pitch. As shown in Fig. 5, the 
forces from the pulse analysis obtained using the deform- 
ing mesh agree very well with the forces obtained using 
the rigidly moving mesh. This good agreement between the 
results verifies the three-dimensional deforming mesh capa- 
bility which was implemented in the code. 

Flutter Results 

Flutter characteristics were determined for seven 
freestream Mach numbers, M c c = 0.499, 0.678, 0.900, 
0.960, 0.990, 1.072, and U41 at a = 0° angle of at- 
tack. Each time-marching calculation was restarted from 
the steady-state solution about the rigid wing, and the mo- 
tion of the wing was initiated by specifying a small initial 
velocity for the first two modes. The resulting transients 
were analyzed for their damping and frequency content with 
the modal identification technique which was previously de- 
scribed. The computed flutter dynamic pressure and fre- 
quency were determined by interpolating the specified dy- 
namic pressures and the computed frequencies to the zero 
damping value of the flutter mode. The flutter mode for all 
freestream Mach numbers considered in this study was dom- 
inated by motion in the first bending mode. A summary of 
the computed flutter characteristics in terms of flutter speed 
index ^ 7 ^ 7 = and nondimensional flutter frequency ratio ^ 
is shown in Table 1. For Moo = 0.9, Ref. 3 reports a com- 
puted flutter speed index of 0.353 and a computed flutter 
frequency ratio of 0.42. These results agree very well with 
those shown in Table 1 for = 0.9 which verifies the 
addition of the structural equations of motion to CFL3D 
version 2 . 1 . 




Reduced Frequency, k 

Figure 5 Comparison of generalized aerodynamic 
forces for the rigid pitch and plunge of 
Wing 445.6 at — 0.9 and a = 0°. 

Table 1 Summary of computed 
flutter results for Wing 445.6. 


Moo 

Flutter Speed 
Index 

Flutter 

Frequency 

Ratio 

0.499 

0.439 

0.597 

0.678 

0.417 

0.539 

0.900 

0.352 

0.425 

0.960 

0.275 

0.343 

0.990 

0.310 

0.373 

1.072 

0.466 

0.541 

1.141 

0.660 

0.764 
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The computed flutter characteristics are compared with 
the experimentally measured values of flutter speed index 
and flutter frequency ratio in Fig. 6. The experimental data 
defines a typical transonic flutter “dip” with the bottom 
near M m = 1 -0. At the subsonic freestream Mach numbers 
(Moo = 0.499 and 0.678), the computed flutter speed indexes 
agree well with the experimental values while the computed 
frequency ratios are slightly larger than the experimental val- 
ues. It is interesting to note that at these subsonic freestream 
Mach numbers the computed flutter results are characterized 
by “hard” flutter crossings. In other words, small changes in 
dynamic pressure result in large changes in the damping of 
the flutter mode. At Moo = 0.9 and 0.96, the computed flut- 
ter speed indexes are less than the experimental values and 
the frequency ratios agree well with the experimental val- 
ues. The computed flutter results at these freestream Mach 
numbers are characterized by a “mild” flutter crossing. Al- 
though there was no experimental flutter point determined 
at Moo = 0.99, computational results are included to aid in 
identifying the bottom of the flutter “dip”. The computa- 
tional results at Moo = 0.99 are compared in Fig. 6 to esti- 
mated values of flutter speed and frequency determined from 
the faired curves in Fig. 16 of Ref. 23. These faired curves 
which were based on the experimentally determined flut- 
ter points, the experimentally determined no-flutter track, 
and analytic calculations are considered to be of reason- 
able accuracy. 24 The computational results at Moo = 0.99 
as well as those at M^ = 1.072 and 1.141 indicate a prema- 
ture rise in the computational flutter boundary as compared 
with the experimental boundary. Although the boundary is 
more sensitive to freestream Mach number in this range, the 
computed flutter results are still characterized by a “mild” 
flutter crossing. 


O Experiment (Ref. 23) 
A Estimated (Ref. 23) 
-Q-Computed 


0.6b 


Flutter 

Speed 

Index 


0.2 
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Figure 6 Comparison of Euler flutter predictions 
with experimental data for Wing 445.6. 


Computational results for this configuration obtained 
with linear theory and TSD methods were previously re- 
ported in Ref. 16. In this report, three sets of flutter results 
were presented.: results from a linear theory subsonic kernel- 
function program, results using the linear potential equation 
and modeling the wing aerodynamically as a flat plate, and 
results using the complete (nonlinear) TSD equation and in- 
cluding wing thickness. Flutter results, flutter speed index 
and flutter frequency, from the subsonic kernel-function and 
the potential equation compare well with the experimental 
data over the range of freestream Mach numbers from 0.338 
to 1.141. ( Note that the subsonic kernel-function results are 
limited to the subsonic freestream Mach numbers.) Results 
from the nonlinear TSD equation for the subsonic freestream 
Mach numbers 0.678, 0.901, and 0.96 indicate that the flut- 
ter speed index is decreased by 1%, 5%, and 19%, respec- 
tively, from the experimental results with a similar decrease 
in the flutter frequency. These results are similar to the trend 
shown by the current computations in Fig. 6. Subsequent un- 
published calculations by the authors of Ref. 16 indicate that 
flutter results for the supersonic freestream Mach numbers 


which were obtained with the nonlinear TSD equation are 
highly non-conservative. This trend is also consistent with 
the results shown in Fig. 6. It is counter-intuitive that the 
higher-order methods utilizing the TSD and Euler equations 
should produce flutter results which compare less favorably 
with the experimental results than with the results based on 
linear methods. However, it is important to note that the 
modeling of the flow physics is incomplete even when us- 
ing the TSD and Euler equations. The existence and effect 
of highly nonlinear flow phenomena such as strong shocks 
and viscous boundary layers must be investigated before any 
conclusions can be drawn. 

Steady-state Mach contours of the initial flowfields on 
the upper wing surface are shown in Figs. 7(a)-7(d) to illus- 
trate the basic flow characteristics at the selected freestream 
Mach numbers of 0.96, 0.99, 1.072, and 1.141 where time- 
marching flutter calculations were made. Mach contours for 
Moo = 0.400, 0.678, and 0.900, which are not shown, indi- 
cate a smooth expansion and recompression of the flow from 
the leading edge to the trailing edge with very little variation 
in the spanwise direction and no supercritical flow. Mach 


7 





Figure 8 Time history of the first generalized displacement 
for the Wing 445.6 at Moo — 0.99 and a = 0°. 


contours shown in Fig. 7(a) indicate that at Moo = 0.96 an 
area of supercritical flow has formed on the wing. This 
area of supercritical flow does not terminate with a shock. 
Figure 7(b) shows that at Moo = 0.99 the areas of super- 
critical flow has expanded, and a normal shock has formed 
near the tip of the wing at approximately 25% of the local 
chord. Mach contours for Moo = 1.702 and 1.141, shown 
in Figs. 7(c) and 7(d), respectively, indicate that with fur- 
ther increases in freestream Mach number, the normal shock 
transitions to an oblique shock located further downstream 
on the outboard portion of the wing at approximately 70% of 
the local chord. Also, at Moo = 1.141, the spanwise extent 
of the oblique shock has increased to approximately 50% 
of the outboard portion of the wing. Rapid changes in the 
steady-state flow conditions from M c c = 0.96 to 1.141 are 
indicated by the formation and movement of a shock at the 
tip of the wing. This range of freestream Mach number also 
corresponds to the range of Mach numbers where the com- 
puted flutter boundary rapidly rises (See Fig. 6). Therefore, 
small variations due to errors associated with modeling defi- 
ciencies and computational deviations could be expected to 
have a large effect on the final flutter speed and frequency. 
Modeling deficiencies could be due to the neglect of the vis- 
cous effects. Computationally, these errors might be due to 
a lack of spatial convergence. 

The rapid changes in the steady-state flowfield condi- 
tions from Moo = 0.96 to 1.141 suggest that, for a given 
freestream Mach number in this range, the flow characteris- 
tics can also change rapidly during an aeroelastic transient. 
Figure 8 shows the time history of the first generalized dis- 
placement qi for a time-marching aeroelastic transient at 
Moo = 0.99 and at a freestream dynamic pressure of 1.12 
times the estimated experimental value. Recall that qi cor- 
responds to the first bending mode which is the dominant 
component of the flutter mode. The damping and frequency 
content of this aeroelastic transient indicates that the wing 
is dynamically unstable at this condition. Instantaneous sur- 


face contours of the pressure coefficient C p are shown in 
Figs. 9 and 10 for the times 7\ and T 2 , respectively, which 
are indicated in Fig. 8. Upper and lower surface pressures 
are shown at these points in time with A C p = 0.02. At 7i, 
contours on the upper surface indicate that an upper surface 
shock has disappeared while contours on the lower surface 
indicate that a lower surface shock has strengthened and 
moved slightly downstream. Similarly, at T 2 , the opposite 
has occurred. The shock has weakened on the lower sur- 
face and strengthened on the upper surface. Figures 9 and 
10 therefore show that during the aeroelastic transient, rapid 
changes in surface pressures occur due to the formation and 
disappearance of a normal shock on the tip of the wing. 
Figures 9 and 10 also illustrate an unusual shock behavior 
during the aeroelastic transient in that there is little chord- 
wise movement of the shock as it strengthens and weakens. 
For two-dimensional airfoils, significant shock weakening is 
usually accompanied by large shock motion. 

Conclusions 

Modifications to an existing three-dimensional, im- 
plicit, upwind Euler/Navier-Stokes code (CFL3D Version 
2.1) for the aeroelastic analysis of wings was described. 
These modifications included the incorporation of a deform- 
ing mesh algorithm and the addition of the structural equa- 
tions of motion for their simultaneous time-integration with 
the governing flow equations. Euler results from calcula- 
tions performed for a rigid wing undergoing forced pitching 
and plunging motions were presented to check the deforming 
mesh algorithm. Aeroelastic Euler results for a 45° swept- 
back wing at a freestream Mach number of 0.9 were com- 
pared to those presented in Ref. 3 to check the addition of 
the structural equations of motion. Calculated flutter results 
for the same 45° swept-back wing were compared with the 
experimental data for seven freestream Mach numbers which 
define the flutter boundary over a range of Mach number 
from 0.499 to 1.14. These comparisons showed good agree- 
ment in flutter characteristics for freestream Mach numbers 
below unity. For freestream Mach numbers above unity, the 
computed aeroelastic results predicted a premature rise in the 
flutter boundary as compared with the experimental bound- 
ary. Steady-state Mach contours of the initial flowfields il- 
lustrated rapid changes in the basic flow characteristics from 
Moo =0.96 to 1.141 which is indicated by the formation 
and movement of a shock near the tip of the wing. Instanta- 
neous surface pressure contours during an aeroelastic tran- 
sient at Moo = 0.99 also demonstrated significant changes 
in the flowfield due to the formation and disappearance of 
a normal shock on the tip of the wing which is induced by 
the aeroelastic motions at this freestream Mach number. 
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